/***********************************************************************
Compile       : gcc -O3 -o subnet subnet.c
Usage         : ./subnet N L M inputfile 
                (./subnet 11 12 2 Budding.DAT)
N             : number of the nodes
L             : number of the steps        
M             : number of the nodes in the subnetworks
inputfile     : solution file generating from program design
Output        : distribution of the subnetworks with size=M
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

struct solution {
	int num;
	int *network;
};

struct subnodes {
	int num;
	int *nodes;
};

struct subnetwork {
	int idx;
	int links;
	int isomer;
	double hist;
};

void print_help ()
{
	printf ("Usage         : ./subnet N L M inputfile\n");
	printf ("                (./subnet 11 12 2 Budding.DAT)\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n");
	printf ("M             : number of the nodes in the subnetworks\n");
	printf ("inputfile     : solution file generating from program design\n\n");
	printf ("Output        : distribution of the subnetworks with size=M\n");
}

int network_to_idx (int N, int network[N][N])
{
	int i, j, idx;

	idx = 0;
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (i == j)
				idx = 2 * idx + network[i][j] + 1;
			else
				idx = 3 * idx + network[i][j] + 1;

		}
	}

	return idx;
}

void idx_to_network (int N, int idx, int network[N][N])
{
	int i, j;

	for (i = N - 1; i >= 0; i--)
	{
		for (j = N - 1; j >= 0; j--)
		{
			if (i == j)
			{
				network[i][j] = idx % 2 - 1;
				idx /= 2;
			}
			else
			{
				network[i][j] = idx % 3 - 1;
				idx /= 3;
			}
		}
	}
}

void print_network (int N, int idx, struct subnetwork subnetwork[])
{
	int network[N][N];
	int i, j, flag;

	idx_to_network (N, subnetwork[idx].idx, network);

	printf ("c subnetwork %d %d %lg\n", subnetwork[idx].idx, subnetwork[idx].links, subnetwork[idx].hist);
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (network[i][j] == 1)
				printf ("g ");
			else if (network[i][j] == 0)
				printf ("n ");
			else if (network[i][j] == -1)
				printf ("r ");
		}
		printf ("\n");
	}
}

void generate_permutation (int p, int M, int P[p][M])
{
	int i, j, k, n, s, select[M], flag[M];

	for (i = 0; i < p; i++)
	{
		n = i;
		for (j = M; j > 0; j--)
		{
			select[M - j] = n % j;
			n /= j;
		}

		s = -1;
		for (j = 0; j < M; j++)
		{
			flag[j] = 0;
		}
		for (j = 0; j < M; j++)
		{
			for (k = 0; k < select[j] + 1; )
			{
				s = (s + 1) % M;
				if (flag[s] == 0)
					k ++;
			}
			flag[s] = 1;
			P[i][j] = s;
		}
	}
}

void generate_isomer (int M, int P[M], int network[M][M], int network2[M][M])
{
	int i, j;

	for (i = 0; i < M; i++)
	{
		for (j = 0; j < M; j++)
		{
			network2[i][j] = network[P[i]][P[j]];
		}
	}
}

void find_isomer(int idx, struct subnetwork subnetwork[], int p, int M, int P[p][M])
{
	int i, j, min, links;
	int network[M][M], network2[M][M], index[p];

	idx_to_network (M, idx, network);

	links = 0;
	for (i = 0; i < M; i++)
	{
		for (j = 0; j < M; j++)
		{
			if (network[i][j] != 0)
				links ++;
		}
	}

	for (i = 0; i < p; i++)
	{
		generate_isomer (M, P[i], network, network2);
		index[i] = network_to_idx (M, network2);
		if (i == 0)
			min = index[i];
		else if (index[i] < min)
			min = index[i];
	}

	for (i = 0; i < p; i++)
	{
		subnetwork[index[i]].isomer = min;
		subnetwork[index[i]].links = links;
	}
}

void set_isomer (int M, int m, struct subnetwork subnetwork[m])
{
	int i, p;

	/* p = M! */
	p = 1;
	for (i = 0; i < M; i++)
	{
		p *= (i + 1);
	}

	int P[p][M];

	generate_permutation (p, M, P);

	for (i = 0; i < m; i++)
	{
		subnetwork[i].isomer = -1;
	}

	for (i = 0; i < m; i++)
	{
		if (subnetwork[i].isomer == -1)
			find_isomer (i, subnetwork, p, M, P);
	}
}

void generate_permutation22 (int N, int M, int subnode[M], int n, struct subnodes *subnodes)
{
	int i, start;

	if (n == M)
	{
		for (i = 0; i < M; i++)
			subnodes->nodes[subnodes->num * M + i] = subnode[i];
		subnodes->num ++;

		return;
	}

	if (n == 0)
		start = -1;
	else
		start = subnode[n - 1];
	for (i = start + 1; i < N; i++)
	{
		subnode[n] = i;
		generate_permutation22 (N, M, subnode, n + 1, subnodes);
	}
}

void generate_permutation2 (int N, int M, struct subnodes *subnodes)
{
	int subnode[M];

	generate_permutation22 (N, M, subnode, 0, subnodes);
}

void do_stat2 (int N, int M, struct solution solution[N], int nodes[M], struct subnetwork subnetwork[], double weight, int select[M], int n)
{
	int i, j;

	if (n == M)
	{
		int idx, network[M][M];

		for (i = 0; i < M; i++)
		{
			for (j = 0; j < M; j++)
			{
				network[i][j] = solution[nodes[i]].network[select[i] * N + nodes[j]];
			}
		}

		idx = network_to_idx (M, network);
		subnetwork[subnetwork[idx].isomer].hist += weight;
		return;
	}

	for (i = 0; i < solution[nodes[n]].num; i++)
	{
		select[n] = i;
		do_stat2 (N, M, solution, nodes, subnetwork, weight, select, n + 1);
	}
}

int compare (struct subnetwork *p1, struct subnetwork *p2)
{
	if (p1->links == p2->links)
	{
		if (p1->hist > p2->hist)
			return 1;
		else
			return -1;
	}
	else
		return p1->links - p2->links;
}

void do_stat (int N, int M, int m, struct solution solution[N], struct subnodes subnodes, struct subnetwork subnetwork[m])
{
	int i, j, select[M];
	double D, D2;

	D = 1.0;
	for (i = 0; i < N; i++)
	{
		D *= solution[i].num;
	}

	for (i = 0; i < subnodes.num; i++)
	{
		D2 = D;
		for (j = 0; j < M; j++)
		{
			D2 /= solution[subnodes.nodes[i * M + j]].num;
		}
		do_stat2 (N, M, solution, subnodes.nodes + i * M, subnetwork, D2, select, 0);
	}
}

int main (int argc, char *argv[])
{
	FILE *fp;
	char buf[256];
	int i, j, k;
	int N, L, M;
	
	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	N = atoi (argv[1]); /* Number of nodes */
	L = atoi (argv[2]); /* Number of steps */
	M = atoi (argv[3]); /* Size of subnetwork */

	int m;

	/* m: number of all possible subnetworks*/
	m = 1;
	for (i = 0; i < M; i++)
	{
		for (j = 0; j < M; j++)
		{
			if (i == j)
				m *= 2;
			else
				m *= 3;
		}
	}

	struct subnetwork subnetwork[m];

	for (i = 0; i < m; i++)
	{
		subnetwork[i].idx = i;
		subnetwork[i].hist = 0.0;
	}

	set_isomer (M, m, subnetwork);

	int m2;

	/* m2 = C(M, N) */
	m2 = 1;
	for (i = 0; i < M; i++)
	{
		m2 *= (N - i);
	}
	for (i = 0; i < M; i++)
	{
		m2 /= i + 1;
	}

	struct subnodes subnodes;

	subnodes.num = 0;
	subnodes.nodes = calloc (m2, sizeof (int) * M);

	generate_permutation2 (N, M, &subnodes);

	struct solution solution[N];

	fp = fopen (argv[4], "r");
	while (fgets (buf, 256, fp))
	{
		for (i = 0; i < L; i++)
			fgets (buf, 256, fp);

		for (i = 0; i < N; i++)
		{
			fgets (buf, 256, fp);
			sscanf (buf, "Total Number of Soultion for Node %*d : %d", &solution[i].num);
			solution[i].network = calloc (solution[i].num, sizeof (int) * N);
			for (j = 0; j < solution[i].num; j++)
			{
				fgets (buf, 256, fp);
				for (k = 0; k < N; k++)
				{
					if (buf[2 * k] == 'g')
						solution[i].network[j * N + k] = 1;
					else if (buf[2 * k] == 'n')
						solution[i].network[j * N + k] = 0;
					else if (buf[2 * k] == 'r')
						solution[i].network[j * N + k] = -1;
				}
			}
		}

		do_stat (N, M, m, solution, subnodes, subnetwork);
	}
	fclose (fp);

	qsort (subnetwork, m, sizeof (struct subnetwork), compare);
	for (i = 0; i < m; i++)
	{
		//if (subnetwork[i].hist > 0.001)
		if (subnetwork[i].isomer == subnetwork[i].idx)
		{
			print_network (M, i, subnetwork);
		}
	}
}
